/***************************************************************************
 *   Copyright (C) 2008 by Chris Bornholdt,,,   *
 *   chris.bornholdt@gmail.com   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#ifndef _CONVOLUTION_HPP_
#define _CONVOLUTION_HPP_

#include <cassert>
#include <cmath>
#include <boost/concept_check.hpp>
#include <boost/concept/requires.hpp>

/*	requirements for ConvolutionKernel	
	double must be convertible to argument_type
	argument_type must be convertible to double
	argument_type/double min_domain() const
	argument_type/double max_domain() const	
*/

template<typename Kernel>
struct ConvolutionKernel : public boost::UnaryFunction<Kernel, 
	typename Kernel::result_type, typename Kernel::argument_type>
{
	public:
		typedef typename Kernel::result_type result_type;
		typedef typename Kernel::argument_type argument_type;
	
		BOOST_CONCEPT_ASSERT((boost::LessThanComparable<argument_type>));
		BOOST_CONCEPT_ASSERT((boost::Convertible<double, argument_type>));
	
		BOOST_CONCEPT_USAGE(ConvolutionKernel)
		{
			support = t.min_domain();
			support = t.max_domain();
		}
		
	private:
		Kernel t;
		double support;
};
	
/*	helper class to perform convolution between two object types discretely.
	specialized to allow proper integer speeds, also will center the kernel 
	over the point of convolution, as though the kernel were defined over the bounds
	-(int)(size / 2) <= x < size - (int)(size / 2) instead of 
	0 <= x < size.  Thus all even sized kernels will be adjusted so that the index size/2
	will be the center and thus be moved towards the positive direction during convolution	*/ 

template<typename SourceIterator, typename KernelIterator, bool SmoothEdges = true>
class discrete_convolution
{
	BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<SourceIterator>));
	BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<KernelIterator>));
	BOOST_CONCEPT_ASSERT((boost::Convertible<
		typename boost::RandomAccessIterator<KernelIterator>::value_type, 
		typename boost::RandomAccessIterator<SourceIterator>::value_type>)); 

	public:
		typedef typename boost::RandomAccessIterator<SourceIterator>::value_type result_type;
	
		discrete_convolution(SourceIterator source_begin, SourceIterator source_end, 
			KernelIterator kernel_begin, KernelIterator kernel_end) : source_begin_(source_begin), 
				source_end_(source_end), kernel_begin_(kernel_begin), kernel_end_(kernel_end)
		{	
			kernel_sum_ = std::accumulate(kernel_begin, kernel_end, 
				result_type(0));
		}

		result_type operator()(int x) const
		{
			int size = std::distance(source_begin_, source_end_);
			assert(0 <= x && x < size);
			
			int n = std::distance(kernel_begin_, kernel_end_);
			int dn = n / 2;
			int i = std::max(-dn, -x);
			int i_end = std::min(n - dn, size - x);

			result_type return_value(0);
			if(i == -dn && i == n - dn)
				for(;i < i_end; i++)
					return_value += source_begin_[x + i] * kernel_begin_[(n - dn - 1) - i];
			else
			{
				result_type sum(0);
				for(;i < i_end; i++)
				{
					sum += kernel_begin_[(n - dn - 1) - i];
					return_value += source_begin_[x + i] * kernel_begin_[(n - dn - 1) - i];
				}
				return_value = return_value * (kernel_sum_ / sum);
			}

			return return_value;
		}
					
		result_type operator[](int n) const
		{	return (*this)(n);	}

		template<typename T>
		operator T ()
		{
			BOOST_CONCEPT_ASSERT((boost::Mutable_ForwardContainer<T>));
			
			T return_value(std::distance(source_begin_, source_end_));
			typename T::iterator iter = return_value.begin();
			
			for(int i = 0; i < std::distance(source_begin_, source_end_); i++, iter++)
				*iter = (*this)(i);
			return return_value;
		}

	private:
		SourceIterator source_begin_, source_end_;
		KernelIterator kernel_begin_, kernel_end_;
		result_type kernel_sum_;
};	
	
template<typename SourceIterator, typename KernelIterator>
class discrete_convolution<SourceIterator, KernelIterator, false>
{
	BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<SourceIterator>));
	BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<KernelIterator>));
	BOOST_CONCEPT_ASSERT((boost::Convertible<
		typename boost::RandomAccessIterator<KernelIterator>::value_type, 
		typename boost::RandomAccessIterator<SourceIterator>::value_type>)); 

	public:
		typedef typename boost::RandomAccessIterator<SourceIterator>::value_type result_type;
	
		discrete_convolution(SourceIterator source_begin, SourceIterator source_end, 
			KernelIterator kernel_begin, KernelIterator kernel_end) : source_begin_(source_begin), 
				source_end_(source_end), kernel_begin_(kernel_begin), kernel_end_(kernel_end)
		{	
			kernel_sum_ = std::accumulate(kernel_begin, kernel_end, 
				result_type(0));
		}

		result_type operator()(int x) const
		{
			int size = std::distance(source_begin_, source_end_);
			assert(0 <= x && x < size);
			
			int n = std::distance(kernel_begin_, kernel_end_);
			int dn = n / 2;
			int i = std::max(-dn, -x);
			int i_end = std::min(n - dn, size - x);

			result_type return_value(0);
			for(;i < i_end; i++)
				return_value += source_begin_[x + i] * kernel_begin_[(n - dn - 1) - i];

			return return_value;
		}
					
		result_type operator[](int n) const
		{	return (*this)(n);	}

		template<typename T>
		operator T ()
		{
			BOOST_CONCEPT_ASSERT((boost::Mutable_ForwardContainer<T>));
			
			T return_value(std::distance(source_begin_, source_end_));
			typename T::iterator iter = return_value.begin();
			
			for(int i = 0; i < std::distance(source_begin_, source_end_); i++, iter++)
				*iter = (*this)(i);
			return return_value;
		}

	private:
		SourceIterator source_begin_, source_end_;
		KernelIterator kernel_begin_, kernel_end_;
		result_type kernel_sum_;
};	

template<bool SmoothEdges, typename Source, typename Kernel>
BOOST_CONCEPT_REQUIRES(
	((boost::RandomAccessContainer<Source>))
	((boost::RandomAccessContainer<Kernel>))
	((boost::Convertible<typename boost::RandomAccessContainer<Kernel>::value_type, 
		typename boost::RandomAccessContainer<Source>::value_type>)),
	(discrete_convolution<typename boost::RandomAccessContainer<Source>::const_iterator, 
		typename boost::RandomAccessContainer<Kernel>::const_iterator, SmoothEdges>))	
discretely_convolve(const Source& source, const Kernel& kernel)
{	
	return discrete_convolution<typename Source::const_iterator, typename Kernel::const_iterator, 
		SmoothEdges>(source.begin(), source.end(), kernel.begin(), kernel.end());	
}

template<typename SourceIterator, typename Kernel, bool SmoothEdges = true>
class general_convolution
{
	BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<SourceIterator>));
	BOOST_CONCEPT_ASSERT((ConvolutionKernel<Kernel>));
	BOOST_CONCEPT_ASSERT((boost::Convertible<
		typename ConvolutionKernel<Kernel>::result_type, 
		typename boost::RandomAccessIterator<SourceIterator>::value_type>));

	public:
		typedef typename boost::RandomAccessIterator<SourceIterator>::value_type result_type;
		
		general_convolution(SourceIterator source_begin, SourceIterator source_end, Kernel kernel)
			: source_begin_(source_begin), source_end_(source_end), kernel_(kernel) {}

		result_type operator()(double x) const
		{
			int size = std::distance(source_begin_, source_end_);
			//allow values outside the boundaries to be used, so long as it is within one unit
			assert(-1 < x && static_cast<unsigned int>(x) < size);
			
			int whole_x = static_cast<int>(x);
			double part_x = x - static_cast<double>(whole_x);
			int i = std::max((int)std::floor(part_x - kernel_.max_domain()), -whole_x);
			int i_end = std::min((int)std::floor(part_x - kernel_.min_domain()) + 1, size - whole_x); 
			
			result_type return_value(0);
			double edge_scale = kernel_.edge_scale(i, i_end);
			for(;i < i_end; i++)
			{
				return_value += source_begin_[i + whole_x] * kernel_(static_cast<double>(-i) + part_x);
			}
			
			return_value = return_value * edge_scale;
			return return_value;
		}
		
		template<typename T>
		operator T () const
		{
			BOOST_CONCEPT_ASSERT((boost::Mutable_ForwardContainer<T>));
		
			int size = std::distance(source_begin_, source_end_);
			T return_value(size);
			for(int i = 0; i < size; i++)
				return_value[i] = (*this)(i);
			return return_value;
		}
		
	private:
		SourceIterator source_begin_, source_end_;
		Kernel kernel_;
};

template<typename SourceIterator, typename Kernel>
class general_convolution<SourceIterator, Kernel, false>
{
	BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<SourceIterator>));
	BOOST_CONCEPT_ASSERT((ConvolutionKernel<Kernel>));
	BOOST_CONCEPT_ASSERT((boost::Convertible<
		typename ConvolutionKernel<Kernel>::result_type, 
		typename boost::RandomAccessIterator<SourceIterator>::value_type>));

	public:
		typedef typename boost::RandomAccessIterator<SourceIterator>::value_type result_type;
		
		general_convolution(SourceIterator source_begin, SourceIterator source_end, Kernel kernel)
			: source_begin_(source_begin), source_end_(source_end), kernel_(kernel) {}

		result_type operator()(double x) const
		{
			int size = std::distance(source_begin_, source_end_);
			//allow values outside the boundaries to be used, so long as it is within one unit
			assert(-1 < x && x < size);
			
			int whole_x = static_cast<int>(x);
			double part_x = x - static_cast<double>(whole_x);
			int i = std::max((int)std::floor(part_x - kernel_.max_domain()), -whole_x);
			int i_end = std::min((int)std::floor(part_x - kernel_.min_domain()) + 1, size - whole_x); 
			
			result_type return_value(0);
		
			for(;i < i_end; i++)
			{
				return_value += source_begin_[i + whole_x] * kernel_(static_cast<double>(-i) + part_x);
			}
			
			return return_value;
		}
		
		template<typename T>
		operator T ()
		{
			BOOST_CONCEPT_ASSERT((boost::Mutable_ForwardContainer<T>));
		
			int size = std::distance(source_begin_, source_end_);
			T return_value(size);
			for(int i = 0; i < size; i++)
				return_value[i] = (*this)(i);
			return return_value;
		}
		
	private:
		SourceIterator source_begin_, source_end_;
		Kernel kernel_;
};

template<bool SmoothEdges, typename Source, typename Kernel>
BOOST_CONCEPT_REQUIRES(
	((boost::RandomAccessContainer<Source>))
	((ConvolutionKernel<Kernel>))
	((boost::Convertible<typename ConvolutionKernel<Kernel>::result_type, 
		typename boost::RandomAccessContainer<Source>::value_type>)),
	(general_convolution<typename boost::RandomAccessContainer<Source>::const_iterator, 
		Kernel, SmoothEdges>))	
convolve(const Source& source, Kernel kernel)
{	
	return general_convolution<typename Source::const_iterator, Kernel, 
		SmoothEdges>(source.begin(), source.end(), kernel);	
}
			
template<typename Source, typename Kernel, bool SmoothEdges = true>
class discrete_area_convolution
{
	public:
		typedef typename Source::value_type result_type;
		
		discrete_area_convolution(const Source& source, const Kernel& kernel)
			: source_(source), kernel_(kernel)
		{	
			kernel_sum_ = std::accumulate(kernel_.begin(), kernel_.end(), 
								result_type(0));	
		} 
			
		result_type operator()(int x, int y)
		{
			assert(0 <= x && x < source_.width());
			assert(0 <= y && y < source_.height());
			
			int n = kernel_.width();
			int m = kernel_.height();
			int dn = n / 2;
			int dm = m / 2;
			int i = std::max(-dn, -x);
			int j = std::max(-dm, -y);
			int i_end = std::min(n - dn, source_.width() - x);
			int j_end = std::min(m - dm, source_.height() - y);
				
			result_type sum(0);
			result_type return_value(0);
			int j_begin = j;
			for(; i < i_end; i++)
				for(j = j_begin; j < j_end; j++)
				{
					sum += kernel_(n - 1 - dn - i, m - 1 - dm - j);
					return_value += source_(x + i, y + j) * 
							kernel_(n - 1 - dn - i, m - 1 - dm - j);
				}
			
			return_value = return_value * (kernel_sum_ / sum);
			
			return return_value;
		}
		
		operator Source()
		{
			Source return_value(source_.width(), source_.height());
			for(int i = 0; i < source_.width(); i++)
				for(int j = 0; j < source_.height(); j++)
					return_value(i, j) = (*this)(i, j);
			return return_value;
		}
		
	private:
		const Source& source_;
		const Kernel& kernel_;
		result_type kernel_sum_;
};

template<typename Source, typename Kernel>
class discrete_area_convolution<Source, Kernel, false>
{
	public:
		typedef typename Source::value_type result_type;
		
		discrete_area_convolution(const Source& source, const Kernel& kernel)
			: source_(source), kernel_(kernel)	{}
			
		result_type operator()(int x, int y)
		{
			assert(0 <= x && x < source_.width());
			assert(0 <= y && y < source_.height());
			
			int n = kernel_.width();
			int m = kernel_.height();
			int dn = n / 2;
			int dm = m / 2;
			int i = std::max(-dn, -x);
			int j = std::max(-dm, -y);
			int i_end = std::min(n - dn, source_.width() - x);
			int j_end = std::min(m - dm, source_.height() - y);
				
			result_type return_value(0);
			int j_begin = j;
			for(; i < i_end; i++)
				for(j = j_begin; j < j_end; j++)
					return_value += source_(x + i, y + j) * 
							kernel_(n - 1 - dn - i, m - 1 - dm - j);
			
			return return_value;
		}
		
		operator Source()
		{
			Source return_value(source_.width(), source_.height());
			for(int i = 0; i < source_.width(); i++)
				for(int j = 0; j < source_.height(); j++)
					return_value(i, j) = (*this)(i, j);
			return return_value;
		}
			
	private:
		const Source& source_;
		const Kernel& kernel_;
};

template<bool SmoothEdges, typename Source, typename Kernel>
discrete_area_convolution<Source, Kernel, SmoothEdges>
discretely_convolve_area(const Source& source, const Kernel& kernel)
{	
	return discrete_area_convolution<Source, Kernel, SmoothEdges>(
			source, kernel);
}

template<typename Source, typename Kernel, bool SmoothEdges = true>
class discrete_separable_area_convolution
{
	public:
		typedef typename Source::value_type result_type;
		
		discrete_separable_area_convolution(const Source& source, 
			const Kernel& kernel) : source_(source), kernel_(kernel)
		{
			kernel_sum_ = std::accumulate(kernel.begin(), kernel.end(), 
				result_type(0));
		}
			
		//TODO replace vector with some stack based solution
		result_type operator()(int x, int y)
		{
			assert(0 <= x && x < source_.width());
			assert(0 <= y && y < source_.height());
			
			int n = kernel_.size();
			int dn = n / 2;
			int i_begin = std::max(-dn, -x);
			int j_begin = std::max(-dn, -y);
			int i_end = std::min(n - dn, source_.width() - x);
			int j_end = std::min(n - dn, source_.height() - y);
			
			std::vector<result_type> horizontal_results(j_end - j_begin);
			result_type result(0);
			result_type sum(0);
			for(int j = j_begin; j < j_end; j++, result = 0, sum = 0)
			{	
				for(int i = i_begin; i < i_end; i++)
				{
					sum += kernel_[n - 1 - dn - i];
					result += source_(x + i, y + j) * kernel_[n - 1 - dn - i];
				}
				result *= (kernel_sum_ / sum);
				horizontal_results[j - j_begin] = result;
			}
			
			result = 0;
			sum = 0;
			for(int j = j_begin; j < j_end; j++)
			{
				sum += kernel_[n - 1 - dn - j];
				result += horizontal_results[j - j_begin] * kernel_[n - 1 - dn - j];
			}
			
			result *= (kernel_sum_ / sum);
			return result;
		}
		
		//TODO rewrite source to convolve whole area in two steps
		operator Source()
		{
			Source return_value(source_.width(), source_.height());
			for(int i = 0; i < source_.width(); i++)
				for(int j = 0; j < source_.height(); j++)
					return_value(i, j) = (*this)(i, j);
			return return_value;
		}
		
	private:
		const Source& source_;
		const Kernel& kernel_;
		result_type kernel_sum_;
};

template<typename Source, typename Kernel>
class discrete_separable_area_convolution<Source, Kernel, false>
{
	public:
		typedef typename Source::value_type result_type;
		
		discrete_separable_area_convolution(const Source& source, 
			const Kernel& kernel) : source_(source), kernel_(kernel) {}
			
		result_type operator()(int x, int y)
		{
			assert(0 <= x && x < source_.width());
			assert(0 <= y && y < source_.height());
			
			int n = kernel_.size();
			int dn = n / 2;
			int i_begin = std::max(-dn, -x);
			int j_begin = std::max(-dn, -y);
			int i_end = std::min(n - dn, source_.width() - x);
			int j_end = std::min(n - dn, source_.height() - y);
			
			std::vector<result_type> horizontal_results(j_end - j_begin);
			result_type result(0);
			for(int j = j_begin; j < j_end; j++, result = 0)
			{	
				for(int i = i_begin; i < i_end; i++)
					result += source_(x + i, y + j) * kernel_[n - 1 - dn - i];
				horizontal_results[j - j_begin] = result;
			}
			
			result = 0;
			for(int j = j_begin; j < j_end; j++)
				result += horizontal_results[j - j_begin] * kernel_[n - 1 - dn - j];
			
			return result;
		}
		
		operator Source()
		{
			Source return_value(source_.width(), source_.height());
			std::vector<result_type> results(source_.width());
			typename Source::iterator iter = return_value.begin();
			
			int n = kernel_.size();
			int dn = n / 2;
			//first, apply convolution in horizontal direction
			for(int y = 0; y < source_.height(); y++)
			{
				int i_begin, i_end;
				for(int x = 0; x < source_.width(); x++)
				{
					i_begin = std::max(-dn, -x);
					i_end = std::min(n - dn, source_.width() - x);
					results[x] = 0;
					
					for(int i = i_begin;  i < i_end; i++)
						results[x] += source_(x + i, y) * kernel_[n - 1 - dn - i];
				}
				
				std::copy(results.begin(), results.end(), iter);
				iter += results.size();
			}
			
			results.resize(source_.height());
				
			for(int x = 0; x < source_.width(); x++)
			{
				int i_begin, i_end;
				for(int y = 0; y < source_.height(); y++)
				{
					i_begin = std::max(-dn, -y);
					i_end = std::min(n - dn, source_.height() - y);
					results[y] = 0;
					
					for(int i = i_begin; i < i_end; i++)
						results[y] += return_value(x, y + i) * kernel_[n - 1 - dn - i];
				}
				
				for(int y = 0; y < source_.height(); y++)
					return_value(x, y) = results[y];
			}
			
			return return_value; 
		}
		
	private:
		const Source& source_;
		const Kernel& kernel_;
};

template<bool SmoothEdges, typename Source, typename Kernel>
discrete_separable_area_convolution<Source, Kernel, SmoothEdges>
discretely_convolve_area_separably(const Source& source, const Kernel& kernel)
{
	return discrete_separable_area_convolution<Source, Kernel, SmoothEdges>(
		source, kernel);
}

template<typename Source, typename Kernel, bool SmoothEdges = true>
class general_area_convolution
{
	public:
		typedef typename Source::value_type result_type;
		
		general_area_convolution(const Source& source, const Kernel& kernel) 
			: source_(source), kernel_(kernel) {}

		result_type operator()(double x, double y)
		{
			assert(-1 < x && x < (double) source_.width());
			assert(-1 < y && y < (double) source_.height());
		
			int whole_x = (int) x;
			int whole_y = (int) y;
			double part_x = x - (double)whole_x;
			double part_y = y - (double)whole_y;
			int i_begin = std::max(floor(part_x - kernel_.max_domain_x()), -whole_x);
			int j_begin = std::max(floor(part_y - kernel_.max_domain_y()), -whole_y);
			int i_end = std::min(floor(part_x - kernel_.min_domain_x()) + 1, source_.width() - whole_x); 
			int j_end = std::min(floor(part_y - kernel_.min_domain_y()) + 1, source_.height() - whole_y);
			
			result_type return_value(0);
			double edge_scale = kernel_.edge_scale(i_begin, j_begin, i_end, j_end);
			for(int j = j_begin; j < j_end; j++)
				for(int i = i_begin; i < i_end; i++)
					return_value += source_(x + i, y + j) * 
						kernel_((double)-i - part_x, (double)-j - part_y);
			
			return_value *= edge_scale;
			return return_value;
		}
		
		operator Source()
		{
			Source return_value(source_.width(), source_.height());
			for(int i = 0; i < source_.width(); i++)
				for(int j = 0; j < source_.height(); j++)
					return_value(i, j) = (*this)(i, j);
			return return_value;
		}
		
	private:
		const Source& source_;
		const Kernel& kernel_;
};

template<typename Source, typename Kernel>
class general_area_convolution<Source, Kernel, false>
{
	public:
		typedef typename Source::value_type result_type;
		
		general_area_convolution(const Source& source, const Kernel& kernel) 
			: source_(source), kernel_(kernel) {}

		result_type operator()(double x, double y)
		{
			assert(-1 < x && x < (double) source_.width());
			assert(-1 < y && y < (double) source_.height());
		
			int whole_x = (int) x;
			int whole_y = (int) y;
			double part_x = x - (double)whole_x;
			double part_y = y - (double)whole_y;
			int i_begin = std::max(floor(part_x - kernel_.max_domain_x()), -whole_x);
			int j_begin = std::max(floor(part_y - kernel_.max_domain_y()), -whole_y);
			int i_end = std::min(floor(part_x - kernel_.min_domain_x()) + 1, source_.width() - whole_x); 
			int j_end = std::min(floor(part_y - kernel_.min_domain_y()) + 1, source_.height() - whole_y);
			
			result_type return_value(0);
			for(int j = j_begin; j < j_end; j++)
				for(int i = i_begin; i < i_end; i++)
					return_value += source_(x + i, y + j) * 
						kernel_((double)-i - part_x, (double)-j - part_y);
			
			return return_value;
		}
		
		operator Source()
		{
			Source return_value(source_.width(), source_.height());
			for(int i = 0; i < source_.width(); i++)
				for(int j = 0; j < source_.height(); j++)
					return_value(i, j) = (*this)(i, j);
			return return_value;
		}
		
	private:
		const Source& source_;
		const Kernel& kernel_;
};

template<bool SmoothEdges, typename Source, typename Kernel>
general_area_convolution<Source, Kernel, SmoothEdges>
convolve_area(const Source& source, const Kernel& kernel)
{	
	return general_area_convolution<Source, Kernel, SmoothEdges>(
		source, kernel);
}

template<typename Source, typename Kernel, bool SmoothEdges = true>
class general_separable_area_convolution
{
	public:
		typedef typename Source::value_type result_type;
		
		general_separable_area_convolution(const Source& source, 
			const Kernel& kernel)	: source_(source), kernel_(kernel) {}
			
		result_type operator()(double x, double y)
		{
			assert(-1 < x && x < (double) source_.width());
			assert(-1 < y && y < (double) source_.height());
		
			int whole_x = (int) x;
			int whole_y = (int) y;
			double part_x = x - (double)whole_x;
			double part_y = y - (double)whole_y;
			int i_begin = std::max(floor(part_x - kernel_.max_domain()), -whole_x);
			int j_begin = std::max(floor(part_y - kernel_.max_domain()), -whole_y);
			int i_end = std::min(floor(part_x - kernel_.min_domain()) + 1, source_.width() - whole_x); 
			int j_end = std::min(floor(part_y - kernel_.min_domain()) + 1, source_.height() - whole_y);
			
			std::vector<result_type> horizontal_results(j_end - j_begin);
			for(int j = j_begin; j < j_end; j++)
				for(int i = i_begin; i < i_end; i++)
					horizontal_results[j - j_begin] += source_(x + i, y + j) *
						kernel_((double)-i - part_x);
						
			result_type result(0);
			for(int j = j_begin; j < j_end; j++)
				result += horizontal_results[j - j_begin] * kernel_((double)-j - part_y);
				
			result *= kernel_.edge_scale(i_begin, i_end) * kernel_.edge_scale(j_begin, j_end);
			return result;
		}
		
		operator Source()
		{
			Source return_value(source_.width(), source_.height());
			for(int i = 0; i < source_.width(); i++)
				for(int j = 0; j < source_.height(); j++)
					return_value(i, j) = (*this)(i, j);
			return return_value;
		}
		
	private:
		const Source& source_;
		const Kernel& kernel_;
};

template<typename Source, typename Kernel>
class general_separable_area_convolution<Source, Kernel, false>
{
	public:
		typedef typename Source::value_type result_type;
		
		general_separable_area_convolution(const Source& source, 
			const Kernel& kernel)	: source_(source), kernel_(kernel) {}
			
		result_type operator()(double x, double y)
		{
			assert(-1 < x && x < (double) source_.width());
			assert(-1 < y && y < (double) source_.height());
		
			int whole_x = (int) x;
			int whole_y = (int) y;
			double part_x = x - (double)whole_x;
			double part_y = y - (double)whole_y;
			int i_begin = std::max((int)floor(part_x - kernel_.max_domain()), -whole_x);
			int j_begin = std::max((int)floor(part_y - kernel_.max_domain()), -whole_y);
			int i_end = std::min((int)floor(part_x - kernel_.min_domain()) + 1, source_.width() - whole_x); 
			int j_end = std::min((int)floor(part_y - kernel_.min_domain()) + 1, source_.height() - whole_y);
			
			std::vector<result_type> horizontal_results(j_end - j_begin);
			result_type tmp(0);
			result_type source(0);
			result_type kernel(0);
			for(int j = j_begin; j < j_end; j++, tmp = 0)
			{
				for(int i = i_begin; i < i_end; i++)
				{
					source = source_(whole_x + i, whole_y + j);
					kernel = kernel_((double)-i + part_x);
					tmp += source * kernel;
				}
				horizontal_results[j - j_begin] = tmp;
			}
						
			result_type result(0);
			for(int j = j_begin; j < j_end; j++)
				result += horizontal_results[j - j_begin] * kernel_((double)-j + part_y);
				
			return result;
		}
		
		operator Source()
		{
			Source return_value(source_.width(), source_.height());
			for(int i = 0; i < source_.width(); i++)
				for(int j = 0; j < source_.height(); j++)
					return_value(i, j) = (*this)(i, j);
			return return_value;
		}
		
	private:
		const Source& source_;
		const Kernel& kernel_;
};

template<bool SmoothEdges, typename Source, typename Kernel>
general_separable_area_convolution<Source, Kernel, SmoothEdges>
convolve_area_separably(const Source& source, const Kernel& kernel)
{
	return general_separable_area_convolution<Source, Kernel, SmoothEdges>(
		source, kernel);
}

#endif
